Main Models and Diagnostics Plots
Arthur M. Albuquerque
July 14, 2021
In this document, I will show our main models and their respective diagnostics plots. The full explanation about these analyses can be found in our preregistered analysis plan: https://osf.io/kvuqr/
In case you are interested in the R code, you can click in the “Code” → “Show All Code” option in the right upper corner of this document.
pacman::p_load(tidyverse, # Data wrangling
brms, # to fit Bayesian models
tidybayes, # to wrangle and plot brms data
rio, # to import/export files
here, # reproducible file paths
metafor) # to calculate log OR
# Load function to plot diagnostics
source(here("final_analyses/functions/diag_plot.R"))
set.seed(123)# Load original data file
d = import(here("final_analyses", "data", "mortality_data.xlsx"))
# Remove studies that have 0 total events in both treatment arms (3 studies)
d =
d %>%
filter((control_events + trt_events) != 0)
# Calculate log odds ratio
d_logOR =
escalc(
measure = "OR", # log odds ratio,
# Tocilizumab/Sarilumab
ai = trt_events,
n1i = trt_total,
# Control
ci = control_events,
n2i = control_total,
data = d
) %>%
as_tibble() %>%
# Dummy approach, tocilizumab = 0, sarilumab = 1
mutate(treatment = ifelse(treatment == "tocilizumab", 0, 1)) %>%
# Calculate standard error
mutate(sei = sqrt(vi)) Sarilumab only model
Model
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu_{sarilumab}, \tau^2)\\ \\ \mu_{sarilumab} & \sim \operatorname{Normal}(0, 1.5^2) \tag{Priors}\\ \tau & \sim \operatorname{Half-Normal}(0.5) \\ \end{align*} \]
## Formula
mf =
# https://bookdown.org/content/3890/horoscopes-insights.html#use-the-0-intercept-syntax
formula(yi | se(sei) ~ 0 + Intercept + (1 | study))
## Priors
priors =
prior(normal(0, 1.5), class = "b", coef = "Intercept") +
prior(normal(0, 0.5), class = "sd") # Half-Normal(0.5)
## Fit model (suppressed)
# model_sarilumab_only =
# brm(
# data = d_logOR %>% filter(treatment == "1"), # only sarilumab
# family = gaussian,
#
# formula = mf,
# prior = priors,
# sample_prior = TRUE,
#
#
# backend = "cmdstanr", # faster
# cores = parallel::detectCores(),
# chains = 4,
# warmup = 2000,
# iter = 4000,
# control = list(adapt_delta = .95),
# seed = 123,
#
# file = here("final_analyses/output/fits/main/model_sarilumab_only")
# )
# Load model
model_sarilumab_only =
readRDS(here("final_analyses/output/fits/main/model_sarilumab_only.rds"))print(model_sarilumab_only)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: yi | se(sei) ~ 0 + Intercept + (1 | study)
## Data: d_logOR %>% filter(treatment == "1") (Number of observations: 9)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~study (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.23 0.17 0.01 0.65 1.00 2875 2635
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.02 0.18 -0.38 0.33 1.00 4396 3463
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.00 0.00 0.00 0.00 1.00 8000 8000
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Diagnostics plots
diag_plot(model = model_sarilumab_only,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
Posterior predictive check
pp_check(model_sarilumab_only, nsamples = 50)
Sarilumab vs. Tocilizumab model
Model
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu &= \beta_0 + \beta_1 x\\ \\ \beta_0 & \sim \operatorname{Normal}(0, 1.5^2) \tag{Priors} \\ \beta_1 & \sim \operatorname{Normal}(0, 1^2) \\ \tau & \sim \operatorname{Half-Normal}(0.5) \\ \end{align*} \]
# Formula
mf =
formula(yi | se(sei) ~ 0 + Intercept + treatment + (1 | study))
priors =
prior(normal(0, 1.5), class = "b", coef = "Intercept") +
prior(normal(0, 1), class = "b", coef = "treatment") +
prior(normal(0, 0.5), class = "sd")
# model_sarilumab_vs_tocilizumab =
# brm(
# data = d_logOR,
# family = gaussian,
#
# formula = mf,
# prior = priors,
# sample_prior = TRUE,
#
#
# backend = "cmdstanr", # faster
# cores = parallel::detectCores(),
# chains = 4,
# warmup = 2000,
# iter = 4000,
# control = list(adapt_delta = .95),
# seed = 123,
#
# file = here("final_analyses/output/fits/main/model_sarilumab_vs_tocilizumab")
# )
model_sarilumab_vs_tocilizumab =
readRDS(here("final_analyses/output/fits/main/model_sarilumab_vs_tocilizumab.rds"))print(model_sarilumab_vs_tocilizumab)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: yi | se(sei) ~ 0 + Intercept + treatment + (1 | study)
## Data: d_logOR (Number of observations: 25)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~study (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.13 0.09 0.01 0.35 1.00 2467 3454
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.18 0.09 -0.34 0.00 1.00 3701 3585
## treatment 0.23 0.16 -0.11 0.52 1.00 5149 4182
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.00 0.00 0.00 0.00 1.00 8000 8000
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Diagnostics plots
diag_plot(model = model_sarilumab_vs_tocilizumab,
pars_list = c("b_Intercept", "b_treatment", "sd_study__Intercept"),
ncol_trace = 4)
Posterior predictive check
pp_check(model_sarilumab_vs_tocilizumab, nsamples = 50)
Sarilumab + Tocilizumab as prior model
Tocilizumab only model
\[ \begin{equation*} \sigma_{i[W]}^2 = \frac{1}{Wa_i+\frac{1}{2}} + \frac{1}{Wb_i+\frac{1}{2}} + \frac{1}{Wc_i+\frac{1}{2}} + \frac{1}{Wd_i+\frac{1}{2}} \end{equation*} \]
where \(a_i\), \(b_i\), \(c_i\) and \(d_i\) represent the number of events (death and no death) in a tocilizumab study \(i\).
W_fun = function(x) {
sew = sqrt((1/(W*x$a + 1/2) + 1/(W*x$b + 1/2) + 1/(W*x$c + 1/2) + 1/(W*x$d + 1/2)))
sew
}# Create data frame with only relevant data
toci_data = d_logOR %>%
filter(treatment == 0) %>%
mutate(a = trt_events,
b = control_events,
c = trt_total - trt_events,
d = control_total - control_events) %>%
select(study, yi, sei, a:d) %>%
as.data.frame()
# Vector of weights, ranging from 0.01 to 1
W = c(0.01, 0.05, seq(from = 0.1, to = 1, by = 0.1))
# Apply W_fun to create a data frame where each new column is the standard error
# per W
weight_data = plyr::adply(.data = toci_data,
.margins = 1,
.fun = function(x) W_fun(x)) %>%
# Select revelant columns
select(study, yi, V1:last_col()) %>%
# Change new columns names
setNames(c(colnames(toci_data)[1:2], paste0("seW_", W)))Now, let’s fit several frequentist random-effect meta-analysis (one per weight)
\[ \begin{align*} y_i & \sim Normal(\theta_{i[W]}, \sigma_{i[W]}^2) \tag{Likelihood}\\ \theta_{i[W]} & \sim Normal(\mu_{tocilizumab[W]}, \tau_{[W]}^2)\\ \end{align*} \]
# Create vector with all standard error column names
column_names =
colnames(weight_data) %>%
data.frame() %>%
rename(weights = ".") %>%
slice(3:n())
# Create empty list to fill below with all fits in the loop
tocilizumab_freq_models_w = list()
# Create empty data frame to fill with overall mean and sd from each fit per
# weight
dat = data.frame(weight = c(0.01, 0.05, seq(from = 0.1, to = 1, by = 0.1)),
mean = NA, sd = NA)
# Run loop to fit all 12 fits and save info in data frame and list
for (i in 1:nrow(column_names)) {
# Pick a column name
weight_value = column_names %>% slice(i) %>% pull()
# Transform data object into the long format to be able to use filter() below
d_w =
weight_data %>%
pivot_longer(seW_0.01:last_col(),
names_to = "sei")
# Filter only data regarding the weight_value above
d_rma =
d_w %>%
filter(sei == weight_value)
# Fit meta-analysis
fit = rma(yi = yi, sei = value,
data = d_rma,
method = "REML", slab = study)
# Save all fits in this list
tocilizumab_freq_models_w[[i]] = fit
# Save overall mean ($beta) and sd ($se) per model in this data frame
dat[i, "mean"] = tocilizumab_freq_models_w[[i]]$beta
dat[i, "sd"] = tocilizumab_freq_models_w[[i]]$se
}
## Save output
# saveRDS(tocilizumab_freq_models_w,
# here("final_analyses/output/fits/main/frequentist_models_prior_fits.rds"))
#
# saveRDS(dat,
# here("final_analyses/output/fits/main/frequentist_models_prior_mean_sd.rds"))Sarilumab (likelihood) + Tocilizumab (prior) model
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu_{sarilumab[W]}, \tau^2)\\ \\ \mu_{sarilumab{[W]}} & = \mu_{tocilizumab[W]} \tag{Priors} \\ \tau & \sim \operatorname{Half-Normal}(0.5) \end{align*} \]
# Same prior for tau in all models
tau_prior = prior_string("normal(0, 0.5)", class = "sd")
# Different priors for mu_sarilumab. 1 per W
# Define priors with string from arguments saved in dat (mean and sd)
# Weight = 0.01
## Create string to input it in prior_string() below
str_001 = paste0("normal(",
dat[1,2], # Mean
",",
dat[1,3], # SD
")")
## Create prior based on string above
mu_prior_001 = prior_string(str_001, class = "b", coef = "Intercept")
# Weight = 0.05
str_005 = paste0("normal(",dat[2,2],",",dat[2,3],")")
mu_prior_005 = prior_string(str_005, class = "b", coef = "Intercept")
# Weight = 0.10
str_010 = paste0("normal(",dat[3,2],",",dat[3,3],")")
mu_prior_010 = prior_string(str_010, class = "b", coef = "Intercept")
# Weight = 0.20
str_020 = paste0("normal(",dat[4,2],",",dat[4,3],")")
mu_prior_020 = prior_string(str_020, class = "b", coef = "Intercept")
# Weight = 0.30
str_030 = paste0("normal(",dat[5,2],",",dat[5,3],")")
mu_prior_030 = prior_string(str_030, class = "b", coef = "Intercept")
# Weight = 0.40
str_040 = paste0("normal(",dat[6,2],",",dat[6,3],")")
mu_prior_040 = prior_string(str_040, class = "b", coef = "Intercept")
# Weight = 0.50
str_050 = paste0("normal(",dat[7,2],",",dat[7,3],")")
mu_prior_050 = prior_string(str_050, class = "b", coef = "Intercept")
# Weight = 0.60
str_060 = paste0("normal(",dat[8,2],",",dat[8,3],")")
mu_prior_060 = prior_string(str_060, class = "b", coef = "Intercept")
# Weight = 0.70
str_070 = paste0("normal(",dat[9,2],",",dat[9,3],")")
mu_prior_070 = prior_string(str_070, class = "b", coef = "Intercept")
# Weight = 0.80
str_080 = paste0("normal(",dat[10,2],",",dat[10,3],")")
mu_prior_080 = prior_string(str_080, class = "b", coef = "Intercept")
# Weight = 0.90
str_090 = paste0("normal(",dat[11,2],",",dat[11,3],")")
mu_prior_090 = prior_string(str_090, class = "b", coef = "Intercept")
# Weight = 1
str_100 = paste0("normal(",dat[12,2],",",dat[12,3],")")
mu_prior_100 = prior_string(str_100, class = "b", coef = "Intercept")
# Put all priors together
priors = list(c(mu_prior_001, tau_prior), # W = 0.01
c(mu_prior_005, tau_prior), # W = 0.05
c(mu_prior_010, tau_prior), # W = 0.10
c(mu_prior_020, tau_prior), # W = 0.20
c(mu_prior_030, tau_prior), # W = 0.30
c(mu_prior_040, tau_prior), # W = 0.40
c(mu_prior_050, tau_prior), # W = 0.50
c(mu_prior_060, tau_prior), # W = 0.60
c(mu_prior_070, tau_prior), # W = 0.70
c(mu_prior_080, tau_prior), # W = 0.80
c(mu_prior_090, tau_prior), # W = 0.90
c(mu_prior_100, tau_prior)) # W = 1.00
# Only data from sarilumab
sari_data = d_logOR %>% filter(treatment == 1)
# Define formula
mf =
formula(yi | se(sei) ~ 0 + Intercept + (1 | study))
# Define arguments for brm()
base_args = list(data = sari_data,
family = gaussian,
formula = mf,
iter = 4000, warmup = 2000, chains = 4,
cores = parallel::detectCores(),
control = list(adapt_delta = .99),
sample_prior = TRUE, seed = 123,
backend = "cmdstanr")
# Create vector with all weights to use later
vector_weights = data.frame(weight = c(0.01, 0.05, seq(from = 0.1, to = 1, by = 0.1)))
# Create list to store all fits
models_sarilumab_w = list()
## Run loop
# for (i in 1:nrow(vector_weights)) {
#
# # Show what model is running
# print(vector_weights[i, c("weight")])
#
# # Put respective prior together with other arguments
# args = c(list(prior = priors[[i]]), base_args)
#
# # Fit model using brm()
# fit = do.call(brm, args)
#
# # Store the whole model in a list, where
# # [[1]] = 0.01 weight, [[2]] = 0.05, [[3]] = 0.10, [[4]] = 0.20, ...
# # [[12]] = 1
# models_sarilumab_w[[i]] = fit
# }
## Change names of list indexing to facilitate data wrangling later
# names(models_sarilumab_w) = vector_weights$weight
## Save list with all models
# saveRDS(models_sarilumab_w,
# here("final_analyses/output/fits/main/models_sarilumab_w.rds"))
# Load models
models_sarilumab_w = readRDS(
here("final_analyses/output/fits/main/models_sarilumab_w.rds"))Priors
vector_weights = data.frame(weight = c(0.01, 0.05, seq(from = 0.1, to = 1, by = 0.1)))
for (i in 1:nrow(column_names)) {
lab = paste0("W = ", vector_weights %>% slice(i) %>% pull() )
forest(tocilizumab_freq_models_w[[i]], transf = exp, xlab = "Odds Ratio")
grid::grid.text(lab, .65, .7, gp=grid::gpar(cex=2))
}
These are the overall estimates from the frequentist models
dat %>%
ggplot(aes(y = "", dist = distributional::dist_normal(mean, sd),
color = sd)) +
stat_dist_slab(fill = NA) +
coord_cartesian(expand = FALSE) +
labs(x = "log odds ratio", y = "Density\n") +
theme_minimal()
dat %>%
ggplot(aes(y = sd, dist = distributional::dist_normal(mean, sd))) +
stat_dist_halfeye(position = "dodge",
.width = .95) +
labs(x = "log odds ratio", y = "Standard deviation of each distribution\n") +
theme_minimal()
Diagnostics plots
Sarilumab + Tocilizumab model
# Run loop
for (i in 1:nrow(vector_weights)) {
# Run custom function diag_plot, 1 per weight
p = diag_plot(model = models_sarilumab_w[[i]],
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
# Display plot
print(p)
# Add legend to show respective weight
lab = paste0("W = ", vector_weights %>% slice(i) %>% pull() )
grid::grid.text(lab, .75, .05, gp=grid::gpar(cex=2))
}
Posterior predictive check
# Run lopp
for (i in 1:nrow(vector_weights)) {
p = pp_check(models_sarilumab_w[[i]], nsamples = 50)
print(p)
lab = paste0("W = ", vector_weights %>% slice(i) %>% pull() )
grid::grid.text(lab, .8, .9, gp=grid::gpar(cex=2))
}
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bayesplot_1.8.0 cowplot_1.1.1 metafor_3.0-2 Matrix_1.3-4
## [5] here_1.0.1 rio_0.5.26 tidybayes_2.3.1 brms_2.15.0
## [9] Rcpp_1.0.7 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7
## [13] purrr_0.3.4 readr_1.4.0 tidyr_1.1.3 tibble_3.1.2
## [17] ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 plyr_1.8.6
## [4] igraph_1.2.6 svUnit_1.0.6 splines_4.0.4
## [7] crosstalk_1.1.1 TH.data_1.0-10 rstantools_2.1.1
## [10] inline_0.3.19 digest_0.6.27 htmltools_0.5.1.1
## [13] rsconnect_0.8.18 fansi_0.5.0 magrittr_2.0.1
## [16] memoise_2.0.0 openxlsx_4.2.3 modelr_0.1.8
## [19] RcppParallel_5.1.4 matrixStats_0.59.0 xts_0.12.1
## [22] sandwich_3.0-1 prettyunits_1.1.1 colorspace_2.0-2
## [25] rvest_1.0.0 ggdist_2.4.0 haven_2.4.1
## [28] xfun_0.23 callr_3.7.0 crayon_1.4.1
## [31] jsonlite_1.7.2 lme4_1.1-27 survival_3.2-11
## [34] zoo_1.8-9 glue_1.4.2 gtable_0.3.0
## [37] emmeans_1.6.1 V8_3.4.2 distributional_0.2.2
## [40] pkgbuild_1.2.0 rstan_2.21.2 abind_1.4-5
## [43] scales_1.1.1 mvtnorm_1.1-1 DBI_1.1.1
## [46] miniUI_0.1.1.1 xtable_1.8-4 foreign_0.8-81
## [49] stats4_4.0.4 StanHeaders_2.21.0-7 DT_0.18
## [52] htmlwidgets_1.5.3 httr_1.4.2 threejs_0.3.3
## [55] arrayhelpers_1.1-0 ellipsis_0.3.2 farver_2.1.0
## [58] pkgconfig_2.0.3 loo_2.4.1 sass_0.4.0
## [61] dbplyr_2.1.1 utf8_1.2.1 labeling_0.4.2
## [64] tidyselect_1.1.1 rlang_0.4.11 reshape2_1.4.4
## [67] later_1.2.0 munsell_0.5.0 cellranger_1.1.0
## [70] tools_4.0.4 cachem_1.0.5 cli_3.0.0
## [73] generics_0.1.0 pacman_0.5.1 mathjaxr_1.4-0
## [76] broom_0.7.6 ggridges_0.5.3 evaluate_0.14
## [79] fastmap_1.1.0 yaml_2.2.1 processx_3.5.2
## [82] knitr_1.33 fs_1.5.0 zip_2.2.0
## [85] nlme_3.1-152 mime_0.10 projpred_2.0.2
## [88] xml2_1.3.2 compiler_4.0.4 shinythemes_1.2.0
## [91] rstudioapi_0.13 curl_4.3.2 gamm4_0.2-6
## [94] reprex_2.0.0 bslib_0.2.5.1 stringi_1.6.2
## [97] highr_0.9 ps_1.6.0 Brobdingnag_1.2-6
## [100] lattice_0.20-44 nloptr_1.2.2.2 markdown_1.1
## [103] shinyjs_2.0.0 vctrs_0.3.8 pillar_1.6.1
## [106] lifecycle_1.0.0 jquerylib_0.1.4 bridgesampling_1.1-2
## [109] estimability_1.3 data.table_1.14.0 httpuv_1.6.1
## [112] R6_2.5.0 promises_1.2.0.1 gridExtra_2.3
## [115] codetools_0.2-18 boot_1.3-28 colourpicker_1.1.0
## [118] MASS_7.3-54 gtools_3.8.2 assertthat_0.2.1
## [121] rprojroot_2.0.2 withr_2.4.2 rdocsyntax_0.4.1.9000
## [124] shinystan_2.5.0 multcomp_1.4-17 mgcv_1.8-36
## [127] parallel_4.0.4 hms_1.1.0 grid_4.0.4
## [130] coda_0.19-4 minqa_1.2.4 cmdstanr_0.4.0
## [133] rmarkdown_2.8 shiny_1.6.0 lubridate_1.7.10
## [136] base64enc_0.1-3 dygraphs_1.1.1.6